{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "In the previous notebook, we showed how QAOA can approximate the ground state $|\\psi_0\\rangle$ of a many-body system characterized by a Hamiltonian $H$. We connected this problem to binary optimization in computer science in notebook 4, and used this connection to understand adiabatic quantum computing and variational algorithms.\n",
    "\n",
    "However, we also talked about the engineering constraints we face in real devices: the evolution in these devices is actually that of an open quantum system, where the quantum processing unit interacts with the environment. In that case, the ground state of $H$ won't be a pure state $|\\psi_0\\rangle$ but a density matrix $\\rho_0$\n",
    "\n",
    "<img src=\"figures/open_system.svg\" alt=\"A quantum processor as an open quantum system\" style=\"width: 400px;\"/>\n",
    "\n",
    "The environment is defined by a temperature $T$, and if we let the system equilibrate, the QPU will become thermalized at temperature $T$. As we saw in the notebook on evolution in open and closed systems, the energy of the states will follow a Boltzmann distribution: $\\rho_0=\\frac{1}{Z} e^{-H/T}$ where $Z=tr (e^{-\\beta H})$ is a normalization factor (called *partition function*), ensuring that $tr(\\rho_0)=1$. If $H$ has a discrete basis of orthonormal eigenstates $\\{|n\\rangle\\}$ with eigenvalues $\\{E_n\\}$, we can write $H=\\sum_n E_n |n\\rangle \\langle n|$ and $\\rho_0=\\frac{1}{Z} \\sum_n e^{-E_n/T} |n\\rangle \\langle n|$ (since exponentiating a diagonal operator consists in exponentiating the elements of the diagonal). Hence, the thermal density matrix is a mixed state where each eigenstate of $H$ with energy $E$ has a classical probability $P(E)=\\frac{1}{Z} e^{-E/T}$, a Boltzmann distribution. We can see that the minimum energy eigenstate will have the highest probability. When $T \\rightarrow 0$, the minimum energy eigenstate will have a probability close to $1$. When $T \\rightarrow \\infty$, all the eigenstates tend to have equal probability.\n",
    "\n",
    "The question that arises now is: how to approximate this thermalized state $\\rho_0$ of the Hamiltonian $H$ using a quantum circuit? For pure ground states, there were two methods: quantum annealing and QAOA. We will see here that those two methods can be adjusted to also prepare thermalized density matrices.\n",
    "\n",
    "We will see later that we can use this preparation to train certain machine learning models.\n",
    "\n",
    "# Quantum annealing\n",
    "\n",
    "Let us start by importing a handful of packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:10:00.429323Z",
     "start_time": "2018-11-19T20:10:00.423825Z"
    }
   },
   "outputs": [],
   "source": [
    "import itertools\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import dimod\n",
    "%matplotlib inline\n",
    "np.set_printoptions(precision=3, suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are interested in the thermal state of the classical Ising model. We create a random model over ten spins and we will sample a hundred states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:10:00.449815Z",
     "start_time": "2018-11-19T20:10:00.437909Z"
    }
   },
   "outputs": [],
   "source": [
    "n_spins = 10\n",
    "n_samples = 1000\n",
    "h = {v: np.random.uniform(-2, 2) for v in range(n_spins)}\n",
    "J = {}\n",
    "for u, v in itertools.combinations(h, 2):\n",
    "    if np.random.random() < .05:\n",
    "        J[(u, v)] = np.random.uniform(-1, 1)\n",
    "model = dimod.BinaryQuadraticModel(h, J, 0.0, dimod.SPIN)\n",
    "sampler = dimod.SimulatedAnnealingSampler()        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's sample the energies at different temperatures. The `dimod` implementation of simulated annealing allows us to set an initial and final temperature for the annealing. If we set it to the same value, we mimic the effect of a finite temperature and we will have a wider range of configurations and energy levels in the samples. The next cell can take a while to execute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:11:12.304140Z",
     "start_time": "2018-11-19T20:10:00.452460Z"
    }
   },
   "outputs": [],
   "source": [
    "temperature_0 = 1\n",
    "response = sampler.sample(model, beta_range=[1/temperature_0, 1/temperature_0], num_reads=n_samples)\n",
    "energies_0 = [solution.energy for solution in response.data()]\n",
    "temperature_1 = 10\n",
    "response = sampler.sample(model, beta_range=[1/temperature_1, 1/temperature_1], num_reads=n_samples)\n",
    "energies_1 = [solution.energy for solution in response.data()]\n",
    "temperature_2 = 100\n",
    "response = sampler.sample(model, beta_range=[1/temperature_2, 1/temperature_2], num_reads=n_samples)\n",
    "energies_2 = [solution.energy for solution in response.data()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define a function to plot the resulting probability distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:11:12.764017Z",
     "start_time": "2018-11-19T20:11:12.307684Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD5CAYAAAAHtt/AAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VNX9//HXmclM9oVAgECEgJRNBUSsGwLu4oJUcN8LVlsr1dpWbeu3v7a22lr3qlgRixbrvu8KLogLslWQXRAJaxYgG8ls5/fHmS0hJJNkZu7M5PN8POYxc+/cufeEh77nzLlnUVprhBBCJA+b1QUQQgjRPhLcQgiRZCS4hRAiyUhwCyFEkpHgFkKIJCPBLYQQSUaCWwghkowEtxBCJBkJbiGESDJpsThpjx49dGlpaSxOLYQQKWvJkiUVWuuito6LSXCXlpayePHiWJxaCCFSllJqcyTHSVOJEEIkGQluIYRIMhLcQgiRZGLSxi2EELHmdrspKyujoaHB6qK0W0ZGBiUlJTgcjg59XoJbCJGUysrKyM3NpbS0FKWU1cWJmNaayspKysrKGDBgQIfOIU0lQoik1NDQQPfu3ZMqtAGUUnTv3r1TvxRiE9x7toCnMSanFkKIgGQL7YDOljs2wV1fAQvujsmphRCiq4tdU8mCu2HHipidXgghuqrYBbfPA6/8DLzumF1CCCGsNGvWLEaNGsWoUaOw2WzB1zfeeGNMrxubXiXK/32w42v47AE4/qaYXEYIIaw0ffp0pk+fztatWzn22GNZvnx5XK4bmxp3bu/Q64/uhPK1MbmMEEIkgpUrV3LYYYfF7XqxCe6cntDncPPa64JXrwOfNyaXEkIIq61YsYJDDz00bteL0QAcBec8BI+OB58byr6CL2fCMdfF5nJCiC6t9JY3Y3bu7+48s81jVq5cySmnnBLc3rhxI3/5y1/Yu3cvL7zwQtTLFLubk70OgXG/Dm3P+zNUfhuzywkhhFWa17gHDhzI448/HrPrxXbk5NgboZf/j/Hsg9d/AT5fTC8phBDx5PP5WL9+PcOGDYvbNWM7V0maE875Jzx2EmgvfLcAljwBR06L6WWFEF1LJM0ZsbJhwwZKSkpwOp1xu2bs5yrpczgc94vQ9nu3QdWmmF9WCCHiYfDgwaxatarJvsrKSq699lqWLVvGHXfcEfVrxmd2wPE3w5o3oGIduOvglZ/ClW+CzR6XywshRDx1796dmTNnxuz88Zkd0JEBP5oJyh/U338Onz0Yl0sLIUSqid+0rn2PgPG/CW1/+BfYsTJulxdCiFQR3/m4j78J+ow2r70ueOknMv2rEEK0U3yD2+6AHz0KaRlme9c38OFf41oEIYRIdvFfAadoMJzyp9D2wvth8+dxL4YQQiQra5YuO/JqGDjBv6Hh5WugscaSogghRLKxJrhtNjOXSXq+2d6zGd79nSVFEUKIZGPdYsH5JXDmP0LbS+fA2rctK44QQiQLa1d5P+w8GD45tP3qdVC93bryCCFEErA2uJWCs+6F3GKzXV8JL/9E5u4WQiSF1Fq6rD2yCuHcf8GcSYCGTZ+YnibH/9LqkgkhRKtSa+my9howrum6lPNvhy1fWVceIYRoh9RYuqwjJtwCJT80r7UXXpwGDXutLZMQQkQgRZYu6wC7A6bMgpljobHadBF840aY8rhpCxdCiAP5f/kxPHfbFcjmS5e98sorvPnmm1RXVzNt2jROPfXUqBYpcWrcAN36w9n3hbZXvgjLn7auPEIIEYHmNe7Jkyfz2GOPMXPmTJ599tmoXy+xghvg0Clw+GWh7bd+DRUbrCuPEEK0orWly26//Xauuy76i6QnTlNJuIl/g++/gMr1ZuGFF66EaR+Yeb2FEKK5CJozYqWlpcu01txyyy1MnDiR0aNHR/2aiVfjBnBmw9TZYPf/Q+xYAe/eam2ZhBCiBS0tXfbggw/ywQcf8MILL8RkJZzErHEDFI+A0/4Kb/3KbC+eDf2Pg8OmWlsuIYRow4wZM5gxY0bMzp+YNe6AI6fDIT8Kbb82A8rXWVceIYRIAIkd3ErB2Q9A4cFm210Hz18BrnpryyWEEBZK7OAGyMiD8+eAPd1s71plmk+0trZcQghhkcQPboDeh8EZd4W2l8+FxY9bVx4hhLBQcgQ3wOjLYeRFoe23b5Ylz4QQXVLyBHdgCtjikWbb54HnLofqbdaWSwgh4ix5ghvAkQkXzIWs7ma7bhc8exl4Gq0tlxBCxFFyBTdAwUFw3r9B2c321sXw5k1ys1II0WUkX3CDmb/71NtD28ueMgN0hBCiC0jO4AY4+qcw4oLQ9ts3m/lNhBAiTrru0mUdpRScdR/sWg07vgaf27R3X/Mx5PWxunRCiC6gay9d1lHOLLhQblYKIawV76XLkrfGHVDQz9ysfHKyWfJs62IzsvLsB2TlHCG6iMPmxC40V1yxou1j4rx0WXLXuAMGjINT/xzaXvqk3KwUQsRN8xr36tWrufbaa5k6dSqPPPJI1K+XGsENcPTP4LDzQ9syslIIESfNa9zDhg1j5syZPPfccyxcuDDq10v+ppIApeDs+6F8tVl4weeGZy+Fq+dBt1KrSyeEiKFImjNi5UBLl7322ms88sgjXHbZZQf4ZMfFpMbt05rdda5YnLp1ziy48OnQzcr6Cnj6AmiwblkjIURqa2npMoBJkybx9ttvM3fu3KhfMybB/c22av7+7ppYnLptBf1MeAeWPStfA89dAV63NeURQqS0lpYu++ijj5gxYwbXXHMNZ5xxRtSvGbOmkh17G2J16rb1OxrOeQheutpsb/wQ3v4NnHmP9DQRQsTchAkTmDBhQszOH7Obk9utDG6AEefD+JtD24tnwxfRv7srhBDxFrPg3lltcXADTLgVDg1bXPjd38Lad6wrjxBCREHMgnt3vZsGtzdWp4+MUqbJpOSH/h0aXvgxbFtmabGEEKIzYtqP29J27gBHhrlZWdDPbLvrYO55ULXR2nIJIUQHxTa4E6G5BCCnCC55ATIKzHZdOfxnCtRVWFsuIUSn6CSdh7+z5U79GndA0RC46BlIyzDbVRtNzdtVZ225hBAdkpGRQWVlZdKFt9aayspKMjIyOnyOmI6cTJgad0D/Y2DKLLNWpfbBtqXw/JX+ft8Oq0snhGiHkpISysrKKC8vt7oo7ZaRkUFJSUmHPx/b4E6kGnfAsLPhjLvMcmcA69+DN26ASf+UPt5CJBGHw8GAAQOsLoYluk5TSbgjp8PxN4W2l/0HPvyrdeURQoh2iGlwb0+0ppJwJ94GIy8ObX/yd5kKVgiRFGIa3DsTtcYNpllk0gMw6OTQvjdvgtWvW1cmIYSIQEyDe1dNAx6vL5aX6By7A86bA30ON9vaZwbobJhnbbmEEKIVMQnuNJu5yefTUFFrwfSu7ZGeAxc/D4UDzbbXBc9cIivGCyESVmyC2x467fa9+2JxiejKKYLLX4U8f/cczz7Tx3tbfFZsFkKI9ohJcDvsoW51CTHZVCQK+pnwzi4y243V8J9zoXytteUSQohmYhTc4TXuJAlugB6D4LJXICPfbNdXwpPnwO7vLC2WEEKEi3lwJ9zoybb0PhQufQkc2Wa7ZjvMmQTV26wtlxBC+MW8qSRhB+G0pmQMXPwM2NPN9p7N8ORkmZRKCJEQYl/jTsbgBhgwDs5/Emz+WQEq1pqat4S3EMJi0lTSmiGnw7n/AuX/e3Z94w/vSmvLJYTo0uLSVJJs0y42cegU+NGjTcP7SQlvIYR1YhLcNqXISTdNDI0eH3vq3bG4TPyMOB8mzwyF986VprdJfZW15RJCdEkxG/LeKy89+Dqpm0sCRl4Akx8B/L8mdq4wNW8JbyFEnMUsuIvzM4Ovk/YGZXMjL2wa3jtWSM1bCBF3Maxxh5blSapBOG0ZdRFMfphQeH8Nc86G2uRbhUMIkZxiWOMOBXdKNJWEG3UxnPMQoWaTlfDERBmkI4SIi4iCWyn1klLqTKVUxEHfKzy4k2GiqfY6/BL4UdgNy8r1MPt0GR4vhIi5SIP4YeBiYL1S6k6l1JC2PlCcF17jbuxY6RLdyAvhvH+Dzb/Q8J7NMHsilK+ztFhCiNQWUXBrrT/QWl8CjAa+Az5QSn2mlLpKKdXi8ui9U73GHTD8HP8q8f5eNDXbTLPJjpXWlksIkbIibvpQSnUHrgSmA8uA+zFB/n5LxzcN7hRr425u8KlwyfOhianqK+DfZ0LZEmvLJYRISZG2cb8MLACygLO11pO01s9qra8Hclr6TGGWE6d/6Ht1g4d6lydKRU5QA8fDZS9Dun9K2IY9pp/3xo+tLZcQIuVEWuN+TGs9XGt9h9Z6O4BSKh1Aaz2mxRPbFD3DB+Gkeq0boN9RcMVrkFlotl21MHcqfPOyteUSQqSUSIP79hb2fd7Wh4q7UnNJQJ9RcNVbkFtstr0ueP4qWPSYteUSQqSMVoNbKdVbKXUEkKmUOlwpNdr/mIBpNmlVr7wU7svdmp7DYNp70P0H/h0a3voVzP8LJPOEW0KIhJDWxvunYW5IlgD3hO2vAX7b1snDa9wpNXoyEgX94MfvwtPnw9bFZt8nf4e6XXDG3WBv659eCCFa1mp6aK3nAHOUUlO01i+29+ThNe6kWTQ4mrK7mzbv566ADf7ON0v+bRZjmDILHJmtflwIIVrSVlPJpf6XpUqpXzZ/tHXy8ImmulyNO8CZDRf9F0ZcGNq35g0zOZWspiOE6IC2bk76OyaTA+S28GhV7/xQr5IuWeMOsDvMrILHXh/at+VLmHUyVKy3rlxCiKTUVlPJo/7nP3bk5L2lxh1is8Gpt0NuH3j3t4CG3ZtMeF/4NJQeZ3UJhRBJotXgVko90Nr7WusZrb3fMzcdpUxHioraRtxeX5P1KLukY35mbly+OB08+/wDdc4xsw2OvMDq0gkhkkBbXRs6NWbbYbfRIyed8ppGtIZdNY30LZAbcgw7y/T1/u+FULsTfG54+SemBj7+ZlCq7XMIIbqsSHqVdErvvAzKa8zsgDv2NkhwB/QdDdM/gLnnQ/lqs++jO6BqE0x6ANLSW/+8EKLLaqtXyX3+59eVUq81f0RygS412VR7FfSDae/CwBNC+75+Bv59FtTstK5cQoiE1lZTyVP+53909AK9u+royUhl5JuZBd/8JSx90uwrWwSPnQAXzoU+h1tbPiFEwmm1xq21XuJ//hgzN8luoAr43L+vTV1mXu7OsDvg7AfgtL+GVtSp3mpW1FnxgrVlE0IknEindT0T+BZ4APgnsEEpNTGSz/buCivhRINScMx1pvYdmBrW0wAvToN5fwKfz9ryCSESRqR98+4GTtBaT9BajwdOAO6N5IPFUuNun0Enw9XzwyaoAhbcDc9cDA17rSuXECJhRBrcNVrrDWHbGzETTbWpV1eeaKqjegwyPU4GnRzat+5t+NcE2PmNZcUSQiSGtnqVnKuUOhdYrJR6Syl1pVLqCuB14KtILhDeVLKruhGfT6Y1jUhmAVz8XNNh8lUbzUjLr5+3rlxCCMu1VeM+2//IAHYC44EJQDkQUYfs7PQ0cjNM5xWX10dVvaujZe16bHYzTH7q7NB6lu56eGk6vPUb8Mi/pRBdUVsDcK6KxkWK8zOoaagFTF/uHjkyuKRdDp0CPYfDs5dBpX9SqkWPwvblcN6/Ia+PpcUTQsRXpL1KMpRS1ymlHlZKzQ48Ir1Ik5VwpJ27Y3oOMzcth50d2rflS3h0HGxaYF25hBBxF+nNyaeA3pgVcT7GrIgT0c1JaNazRAbhdFxGHpz/FJzyp1B/77pys5r8R38Dn9fa8gkh4iLS4B6ktb4NqPPPX3ImcFSkF+ktNe7oUQqO+wVc/ipkF5l92gcf/dXMMlizw9ryCSFiLtLgdvuf9yilDgXygZ6RXiR8Xm6pcUfJgHFwzSfQP2we7+8WwCPHwYYPrCuXECLmIg3ufymlugG3Aa8Bq4C/RXqRYploKjby+sDlr5mpYPFPBVtfAf+ZAu//AbzuVj8uhEhOEQW31nqW1nq31vpjrfVArXXPwOo4keglE03Fjj0NTvitWZQ4p3do/8L74IkzzDSxQoiUEmmvku5KqQeVUkuVUkuUUvcppbpHehGpccfBgHFw7adw8EmhfWWLYObxsPxpswyRECIlRNpU8gywC5gCTAUqgGcjvUhBloNMhx2A2kZPcGEFEWU5RXDJC3DyH8Hm76LvqoFXfgrPXwH1VdaWTwgRFZEGd7HW+s9a603+x+1Ar0gvopRicO/QovBrd0Tck1C0l80GY2+Aae9B90Gh/atehUeOhW/nW1c2IURURBrc7ymlLlRK2fyP84F323OhYWHBvWZHdXs+Kjqi7xGm18mYH4f21WyHp34E79wKbmmyEiJZtTXJVI1Sqhq4GngacPkfzwA/ac+FhoYF9+rtUuOOC2c2nHUvXPQMZPUI7f/iYTPisqxTa0ELISzS1go4uVrrPP+zTWud5n/YtNZ57bnQ0OLQ4VLjjrMhE+Fnn8MPTgvtq1gLj59sug1K7VuIpBJpUwlKqUlKqX/4H2e190LhNe71u2rxeGVFl7jK6QkXPwtn3hOaaVD7TLfBR8dB2WJryyeEiFik3QHvBH6BGXizCviFUuqO9lyoIMsZ7Bbo8vj4rrKunUUVnaYUHDkNfvYZlB4f2l+xFh4/Bd7/P6l9C5EEIq1xnwGcorWerbWeDZyOma+kXYZIO3di6FZqRlyeeXez2vf98OjxUvsWIsFF3FQCFIS9zm/twIp9FTy39jne2fQOC7cu5Ovyr9m0dxP9e/pAeQBp57aczQZHTjdt3wPGhfZXrDO177dvhkb5chUiESkdwYg6pdRFwJ3Ah5hJMcYBt2itWxyEkzkgUw/6f4NaegsA7UvDqbIpKSgk15lrHo7c0GtnLnnOvCbb4fvS7bIQQ1RpDUuegPduA1dtaH9eXzjjLhja7h9XQogOUEot0VqPafO4toJbKaUw8297gCP9uxdprQ84f2hbwd1ZTptzv1DfL+Bb+SJIt6dj/izRxO7N8MYN+w/SGXoWTPw75Pe1plxCdBFRC27/yVZorQ+L9OIHDT9IT5s1jRpXTejhrqG6sZqqhr0oZW2PkjRbWijII6zp5zhygtuZaZmpG/xaw4oX4J1bzEyDAc5cOOk207xis1tXPiFSWLSDew7wT611RCu7jxkzRi9e3PINrtPv/4Q1OypR9gbuvWgIpUU2ql3VTUPeH/TN91W7qql2VePxeSIpRsykqbQD1/TDtnMcOS1+EWSlZSV+8NdXwQd/gKVPNt3fZ7S5qdl3tDXlEiKFRRrcrS4WHOYo4FKl1HdAHaadW2utR7S3YMN757Fmew3a46SmupBRh5S26/Naaxq9jU3CfL+avauaWlftfqEfeHb7OjdPtUd72N24m92Nuzv0eZuyke3IJteRS44zJ1ibb/I6rIYfvi/HmROf8M8qhEkPwogLTfNJxTqzf9tSeOxEGH05nPQHyI54kkghRJREGtyntX1IZJp0CezAZFNKKTLSMshIy6Aoq6hDZQgEf0s1/fB9ta5aqt37H9Po7dzshj7tC56LDnZnP1D4Nwn6Fr4I2h3+pceZ6WIX3g+f3AVeF6Bh6RwzcdWJvzfzoUjziRBx02pwK6UygGuBQcAK4HGtdafaKZoMfd9uTZfAdHs66Znp9Mjs0fbBLQiv8R8w9N21LX4x1Lpr2efZ1+m/IZrhH2jDz3HmHPhXQL/DyJn6ELlL5pC7+QtyfJrcxr1kvvUr1JI5pvdJ/2M6/XcJIdrWVo17Dma9yQXARGA4ZgRlhw1rNr2rz6ex2RK8vbeZzga/2+em1lVrmnPcNU2ea921wdp+8LXbHFvtqg6+bvB2foRjk/Bvj5I+wZd2rcn2VZH7/o/JSc8jp9vB5GZ2bzH885x5Lf4KSOmbvULEQFvBPTzQm0Qp9TiwqLMXLMpNpzDbSVWdizqXl6179nFQYVZnT5tUHDYH3TK60S2jW4fP0Z7wD7T5N/8iiEb4e5Wi2m6n2g746qFyRbvPEaj55zhyTPOPMze4HQj6SLblC0B0FW0Fd/AuntbaE43/KZRSDO2dy2ffVgKwent1lwvuaIhK+HvdwRBvHv6BG71t/QqwtObfjF3ZQ18ATtP+3+K2Myd4fyDb6f8C8H8JZDuyybBnyBeASGhtBfdI/3zcYHqSZPq3A71K2jW1a8DQ3nnB4F6zo4ZTD+ndxidELDjsDrrZoxT+386jZsHfqdnzHbU2GzU2m3nO601t6XHUZhXENPwBvNob7DLa0XZ/MN09A4EeXrNvq9Yf3OffloFeIlZaDW6tdUy6CgwtltVwUkUw/A+ZCsN+BMv+Ax/+BWp3mgOqa6BsPQw+HU75ExQN2e8cbp+bOlcdte5a6tx1wUAPbNe4aoL7W9vubG+fAI/2sLdxL3sb93bqPOFfAC01A7X1KyBwnHwBiOYi7Q4YVeFzc6+RWQJTh80OR1wBh04x3Qc/exACPWjWvQPr34fRl8G43zQZPu+wOSjIKKAgo+AAJ46M2+s2Qe72B7s//GvdtcEvhsCXQvALovm2qxaXz9WpcgRE6wsg0AQU/shx5JDlyAoGfKTvy5dAarAkuH/QMxebAp+GTZV17HN5yXRKP+CUkZ4DJ/4OjrjS1L6XPw1o0F5Y8m9Y/l/44dUw9pdRHcDjsDsosHf+C8Dlde0f/BHU/Jv/SujsQK+AJk1AnZSm0oKB3pHgD+535uC0OeVLwCIRDXlvr9aGvAecePdHbCw3DZGv/fw4RpR07n82kcC2/w/e+z1s+qTpfmcuHHOdeWR06HZJQnN5XcHafou/AiL4VVDnrovaL4BoCzQFZadl7/ec48whK80EfJv7HdnyJeAX7SHvUTesd14wuNdsr5HgTmXFI83CDRs/gnl/MsPmAVw18PGdsOhfMPZGUwt3ZFpa1Ghy2p0U2gspzCjs1HkCTUCBsG/+qHXXUu+uj+j9aP0KgGZNQZ1c0CrNlrZfTb+1Gn92mtnOcmQ12e4q3UItC+6hvXN5c8V2AFbLDcrUpxQcfAIMnABr3oD5t0P5GvPevip4/zaz+vzYX5p5UBwZVpY2oUSrCQhCzUAdDf7wfdGc7M3ji879AACFCoZ5MNgd2cHtbEc2WWlZweAP/2LISstqsi8rLSshvwisC+4mQ9/lBmWXoRQMOxuGnAFfPwcf/RX2fG/eq9kOb/8aFtwNx/3CtJE7pY9/NDntTpx2Z6e6gAYE7wV0IPibvx/NLwGNDp6bzs8ugU3ZTKA3D/WWvgTCav5Zjqzge+HHR2OcgKU17oA1O6rRWifct5qIIZsdRl1keqAsnWMmsAp0IazdAe/eCp/eA8fOMAscO7OtLa/YT7S/BDoS/PXueuo99U22ozUuIMCnfcF7D1H9Igj7JRAI9khZFtwl3TLJSU+jttHD7no3u2oa6ZUnP4+7nDSnadsedYkJ8E/vM8ENUFdumlAW3gfH/Nwcl57b+vlEUorW/QAwzS77PPuCQV7nrqPOExb07vom28HQ99Tvtx3N8QEBTb4IOsiy4FZKMaR3Lks2mzmtV2+vluDuypxZcPRP4YirYNlT8Om9UL3VvFdfCfP+CJ89AEddCz/8iZkvXIgWpNlCC51Eg8fn2S/UW6rth4d9i18E/tfR6CVkWXCDaS4JBPfaHTVMGNLTyuKIRODIMDXr0ZfD8rmw4B7Yu8W8t283fHSHGdxz+GWmG2G3/taWV6S8wFKHec7odFl1+9yhmn+zXwOncEpkZYpKSTqoyQ3KDiyqIFJYWrpZoGHUpfD1M/DJP2DPZvOeux4WPQpfzYJDzzU3MntHvCSqEJZy2Bzkp+eTn57f4XPYoliedgufm3u1RYsqiASX5jS17+uXwLmzoFdYQGsvrHgeZo6Fp86FjR+bxY6FSHGWBvfgsOD+trwWl8fa1d9FArM7YMR5cO0CuPRFGDCu6fvfzoMnJ5kQX/okuKNw+1+IBGVpcOdlOCjpZkbKub2ajRUdv8squgilYNDJcMXrcPWHMHwyqLD/jHeuhNeuh3uGm1Ga1dusK6sQMWJpcIPMFCg6oe9oOH8O/HwxjJkG4f1g91WZgTz3HQYv/Bi2fGVdOYWIsgQI7tANShn6Ljqk+8Fw1j3wy1Vwyp8hv1/oPZ8HVr4Ij58Mj50I/3sW3NEdoCFEvFke3If0CQX3F/5VcYTokMxucNwMmLEMzn8K+h/X9P2tS+Dln8A9w+Dd30Hlt9aUU4hOsjy4jz24B3b/Ku//K9vLrhqpDYlOsqfB8Elw1VtwzSdmVKbdGXp/XxV8/k94cDTMmQTfvALe6M2aJ0SsWR7c+VkOjiwNzXXw4ZpdFpZGpJzikTD5YbhxFZz4+6bNKACbPobnrwjdzNy92ZpyCtEOlgc3wElDewVfz1stwS1iIKcIxv0afrEcLn7ezE4Y3hulbpe5mXn/SHjyHPj6eelSKBJWYgT3sNBQ9wXrK2hwey0sjUhpNjsMPhUu+i/csALG3wy5xWEHaLPgw0vT4R9D4PUbTI8UGdgjEkhCBPfAohwG9DDTdu5ze/l8o9ykFHGQXwIn/BZuWAkXzIWDTwLCphZu3AtLnjA9Uh46ysxcWLPDsuIKEZAQwQ1w0tBQrXu+NJeIeLKnwbCz4LKX4MaVpi28cGDTYyrWwgd/MD1S5p4HK14AVyfX6xKigxInuIeFt3PvJBaLGAvRpvwS0xZ+/VK46m0zyZUjbBEH7YP178GL0+CuQfDidFj3rvRKEXGVMME9prQbuRlmssJtextYLaMohZWUgv7HwuSH4FfrYPIj0H9s02Pc9WaSq6fPh38MhjduhM2fgU/m3BGxlTDB7bDbmszHPX/NTgtLI0SY9BwYdTFc9SbMWG6aUnoMaXrMvipYPBuemGiG2b//f7BtudzUFDGRMMENTdu5P5B2bpGICgeYppTrvoRrPzVzgeeVND2musws9vCv8aZ74Xu3QdliCXERNSoWbcljxozRixcvbvfn9tS7GP3n9/Fp80t10W9Ppig3PerlEyKqfD7Y8oVpNvnmZbNST0vySswK98PPgYOOAltC1ZtEAlBKLdFaj2nruIT6L6cgy8mY/mYtQa3hw7VS6xZJwGYz7eFn3Qs5WUZMAAAMRUlEQVQ3rYOLn4MRF0B6s6Wuqsvgy0fgidPhnqHw5k2mz7jc2BTtZOnSZS05aVhPFn1XBZjeJeePOcjiEgnRDmlOGHyaeXgazao8q16FtW82rYnX7jRLr301C9Lz4Qcnw+CJ5jmz24HPLwQJ1lQCsGFXLSff8zEAWU47y/7vFNLT7NEsnhDx53XDdwtMiK9+A+orWj5O2U3tffDpMGSimbJWdBmRNpUkXHBrrZnwj4/YXFkPwJwf/5Dxg4uiWTwhrOXzmm6Dq1+DtW+HVrFvSY/BoRAvOdIs4SZSVqTBnXBNJUopThrai9kLNwEwf/VOCW6RWmx2GHC8eUz8u1lube07sPYt2La06bEV68zjswdMm/nA8WZo/qCToKBfy+cXKS/hatwACzdUcMmsLwHoW5DJpzefgFKqjU8JkQJqdpiRmGvfNjcuPa3MUNhjsFl/c9BJZtEIR2bciiliI2lr3ABHlhaSm55GTaOHrXv2sXZnTZMlzoRIWbm94YgrzMNVb+YLX/s2bJhneqWEC9TGv3gY0jJMeB98oqmV9zxEuhumsIQMbmeajXGDi3hzxXbAzNEtwS26HGeWadseMtH0j61YBxs+MCG+eSF4wlaL8jTAt/PMAyCz0DTFlB4PA8ZDjx+YwREiJSRkcIPpFhgI7g9W7+S6EwZZXCIhLKQUFA0xj2OuM4s8bF4IG+abMK9Y2/T4fVWmB8uqV812Tm8YMC706NY//n+DiJqEDe4JQ3piU+DTsOz7PazaVs3wPlLrFgIw7dmDTjYP/gp7tpja9qZPzKOuvOnxtTtgxXPmAebGZv+x0P8Y6HcMdB8kNfIkkpA3JwOum7s0WOs+9/C+3HPBqE6fU4iUpzWUrwmF+HcLoGFv65/J6gH9jjZ9yPsdDb1HmnnKRVwlbT/ucP/bsodzHloIQJpNseDmEyjOlzvnQrSLzws7vg4F+ebPzJS0rXFkQ8mYUJCXHAnO7NY/IzotJYIb4PxHP2fRJjME/ppxA7n1jGFROa8QXZbHBduWwfef+x9fQMOe1j+jbNBzuAnzvmPMc48h0nMlylImuD9YtZPpT5pz5aan8dmtJ5KbIaPHhIgan880rYQHeWujOQPS86DP4U3DPKdn258TB5TU/bjDnTi0JwOLstlYXkdNo4dnFm3h6nED2/6gECIyNhv0Gm4eR04z+/ZsMQH+/WfmeddqoFklr7Ha9DPf9HFoX0E/6HsEFI+CPqOg9wjIKozbn9JVJHyNG+C/i77n1pdWAFCcn8EnvzkBh11+ogkRNw3Vpnll62IoWwJlX0FdhNMuF/SH4pEmyItHmlDP7hHb8iaplGkqAWhwexn7t/lU1LoAuO+CUUw+vG/Uzi+EaCetTXNK2WLY6g/y7f9rOiioNXklTcO816GQ16fLd0lMqeAGeHDeeu5+fx0AQ3vn8sb1Y0mTWrcQicPrNhNmbVtmQnzbcti1CryuyD6fUQC9Dgl7HApFQ82an11EygX37joXx945n31uLwC/Pm2IjKYUItF5XFC+OhTk2/9nwj3SmjlAtwH7B3q3UjPLYopJueAGeOjDDdz1rhna67ArXvv5WIYVy2hKIZKK12OG6AeCfMfXsHMVNLYxSCicPd3Mv9JjcGgqgB5DzMITacm7Tm1KBrfH62PqzM9ZvsX0OR1WnMer1x2HM02aTIRIalrD3jLY+Y2pke9aZV5XrAftjfw8yg6FA0yIFw02TS09BptHEjS5pGRwA3xbXssZ9y+g0eMD4PoTB3HTqUNici0hhMXcDWZWxOaBXruz/efK62uaXQr9j25hz5kF0S97B6RscAPM/nQTf3pjFQB2m+LFnx7LqIMS4x9eCBEH+/aYQC9fA+VrQ6/3fN+x82UWhoX5wKbBntMrbr1dUjq4fT7NRY99wZf+ofAHF2Xz5ozjyXCk3s0KIUQ7uOqhcr0J8/K1Jswr1kHVRvB5OnZOR7a5GVo4wP8cFuz5B0V1Mq6UDm6ALVX1nH7fJ9S5TPvXtLEDuO2s4TG9phAiSXlcpja+exNUbQo9V22E3d+Bt7Fj57WlmdGi3Zo1vxQONCHfzuXkUj64oemISqXgv1cfzdEDu8f8ukKIFOLzQc32ZmEeFvBtTYnbmhtWQsFBER+eMnOVtObCIw/inZU7+HhdOVrDr57/H6//fCzdsp1WF00IkSxsNsjvax6lY/d/v76qaS09POBrd7RyXocZDRoDSR3cSin+NmUEp977MdUNHsp27+OSWV8yd/pREt5CiOjIKjSPkiP2f89Vb5paWqqtp6XHbJBQUjeVBLy1YjvXPb2UwJ8yvDhPwlsIYS2fr93zlUfaVJISI1fOOKyYv08ZEeyxs2p7NRfP+pLK2g7ecBBCiM6K4SITKRHcAOeNOYi7po4Mhvfq7dWc9eCnLP1+t7UFE0KIKEuZ4AaYekRJk/DevreBCx79nCcWbiIWTUJCCGGFlApuMOE9+4ojyc80y5u5vZo/vr6K6XMWU7a7jQVShRAiCaRccAOcMLQnb1w/lhEl+cF989bs4pR7PuHBeevZu89tYemEEKJzUjK4AQ4qzOL5a4/hymNLg/v2ub3c/f46jrtzPne8tZqd1e2YE1gIIRJESnQHbMuSzVX87uWVrNlR02S/027j3NF9ueyY/gwvzkN18WWThBDW6hJD3tvD7fXx6vJtPPrxt6zfVbvf+wcXZXP2yD6MH1zEiJIC7DYJcSFEfElwH4DPp5m3ZhcPf7SBZd/vafGYvIw0jj24B0cNLOSQPvkMLc4lL8MR55IKIboaCe42aK1ZtKmKuV9+z/urdgbXsjyQfoVZDC/OY1hxHqU9sjioMIuDumXRI8cpTSxCiKiQ4G6HepeH+Wt28eGacj7dUM7O6shHXKan2eiRk06P3HSKcpzmdU46BVkOctLTyMlIM89hr7OdaaQ7bDjtNlmpXggRJMHdQVpr1u+qZeGGClZurWbV9mrW76zB44vNAB67TeG024JBHni22xQ2pbDbVNPXSmGz0cI+82y3mdeBJvrAbwGlVPA1wfdUcLCSgrDXyv+Z8IU/wvYFz2mOPaJ/NyYf3jcm/z5CdCVdYlrXWFBKMbhXLoN75Qb3NXq8bNhVy6pt1azdUcOW3fV8X7WPsqp6aho7uKqGn9en2efzttlUk8jcXp8EtxBxJMEdgfQ0O4f0yeeQPvlN9mutqW30UFnroqK2kYraRsprXVTUNFLT4KG20U1to4eaBg91jR5qGz3UNphnl9dHo8eHjMQXQrSXBHcnKKXIzXCQm+GgtEd2uz+vtcbj07g8JsTNsxeXx4dXa7w+jc9H6HVwn262j6bva43WoNH+6xD8gtBh1w5+ZzQ7NnBc6HXYeUKFD74eVJTT7r9dCNFxEtwWUkrhsCscdhvZ6VaXRgiRLKRLgxBCJBkJbiGESDIS3EIIkWQkuIUQIslIcAshRJKR4BZCiCQjwS2EEEkmJnOVKKXKgc1RP7EQQqS2/lrrorYOiklwCyGEiB1pKhFCiCQjwS2EEElG5ioRSUUp5QVWhO16Rmt9p1XlEcIK0sYtkopSqlZrHdXpCJVSaVrrzk2sLkQcSVOJSAlKqe+UUn9USi1VSq1QSg31789WSs1WSi1SSi1TSp3j33+lUuo1pdR8YJ5SyqaUelgptUYp9b5S6i2l1FSl1IlKqVfCrnOKUupli/5MIQAJbpF8MpVSy8MeF4S9V6G1Hg08AvzKv+93wHyt9Q+BE4C7lFKBydNHA1O11uOBc4FSYDhwGXCM/5gPgaFKqUAXrauA2TH624SIiLRxi2SzT2s96gDvveR/XoIJYoBTgUlKqUCQZwD9/K/f11pX+V+PBZ7XWvuAHUqpDwG01lop9RRwqVLqCUygXx69P0eI9pPgFqmk0f/sJfTftgKmaK3Xhh+olDoKqIvwvE8ArwMNmHCX9nBhKWkqEanuXeB6pcx69Uqpww9w3EJgir+tuxcwIfCG1nobsA34PSbEhbCU1LhFsslUSi0P235Ha31LK8f/GbgP+FopZQM2AWe1cNyLwEnAKmALsBTYG/b+XKBIa726M4UXIhqkO6AQfkqpHK11rVKqO7AIOE5rvcP/3j+BZVrrxy0tpBBIjVuIcG8opQoAJ/DnsNBegmkPv8nKwgkRIDVuIYRIMnJzUgghkowEtxBCJBkJbiGESDIS3EIIkWQkuIUQIslIcAshRJL5/zRkhcA1RfuQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_probabilities(energy_samples, temperatures):\n",
    "    fig, ax = plt.subplots()\n",
    "    for i, (energies, T) in enumerate(zip(energy_samples, temperatures)):\n",
    "        probabilities = np.exp(-np.array(sorted(energies))/T)\n",
    "        Z = probabilities.sum()\n",
    "        probabilities /= Z\n",
    "        ax.plot(energies, probabilities, linewidth=3, label = \"$T_\" + str(i+1)+\"$\")\n",
    "    minimum_energy = min([min(energies) for energies in energy_samples])\n",
    "    maximum_energy = max([max(energies) for energies in energy_samples])\n",
    "    ax.set_xlim(minimum_energy, maximum_energy)\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    ax.set_xlabel('Energy')\n",
    "    ax.set_ylabel('Probability')\n",
    "    ax.legend()\n",
    "    plt.show()\n",
    "\n",
    "plot_probabilities([energies_0, energies_1, energies_2], \n",
    "                   [temperature_0, temperature_1, temperature_2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just as we saw in the introduction, the distribution flattens out at a high temperature ($T_3$). On the other hand, the energy is peaked for a low temperature, and we do not even have samples for high-energy states."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantum Approximate Thermalization\n",
    "\n",
    "There are many results for preparing a thermal state on a gate-model quantum computer, but most of them need a large-scale device. More recently, a protocol for approximating thermalization was developed using shallow circuits [[1](#1)]. The idea is that if we knew that the thermal state was a pure state $\\psi$ (which means $\\rho=|\\psi \\rangle \\langle \\psi |$), we could apply QAOA to get to the thermal state of a target Hamiltonian. Since QAOA approximates the adiabatic pathway, it should be a conservative change, so at the end of it, we would be close to the thermal state of the target Hamiltonian.\n",
    "\n",
    "To find the thermal state of the simple system, the trick is to purify $\\rho$ on a larger Hilbert space. If we call $\\mathcal{H_1}$ our current Hilbert space, purifying a density matrix $\\rho$ consists of finding a second Hilbert space $\\mathcal{H_2}$ such that there exists $| \\psi \\rangle \\in \\mathcal{H_1} \\otimes \\mathcal{H_2}$ such that $\\rho = \\textrm{Tr}_{\\mathcal{H_2}} \\left( |\\psi \\rangle \\langle \\psi | \\right)$, where $\\textrm{Tr}_{\\mathcal{H_2}}$ is the partial trace taken over the second Hilbert space -- in essence, we are marginalizing the probability distribution. This resembles the idea of what we shown in the notebook on measurements and mixed states: if we trace out a subsystem of the maximally entangled state $|\\phi^+\\rangle$, we get the maximally mixed state. The maximally mixed state is essentially a thermal state at infinite temperature.\n",
    "\n",
    "It can be shown that $| \\psi \\rangle =\\sqrt{2 \\cosh 1/T} \\sum_{z \\in {-1,1}} e^{- z/T} |z \\rangle_{\\mathcal{H_1}} \\otimes | z \\rangle_{\\mathcal{H_2}}$ purifies $\\rho=\\frac{1}{Z}e^{- H_m/T}$ [[1](#1)], where $H_m$ is the mixing Hamiltonian in QAOA. This state can be built with a circuit composed uniquely of RX gates and CNOT gates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:11:13.926133Z",
     "start_time": "2018-11-19T20:11:12.766644Z"
    }
   },
   "outputs": [],
   "source": [
    "import itertools\n",
    "import numpy as np\n",
    "from functools import reduce\n",
    "from qiskit import Aer, BasicAer, QuantumRegister, QuantumCircuit, ClassicalRegister\n",
    "from qiskit import execute\n",
    "from qiskit.quantum_info import Pauli\n",
    "from qiskit.aqua import QuantumInstance\n",
    "from qiskit.aqua.operator import Operator\n",
    "from qiskit.aqua.components.optimizers import COBYLA\n",
    "from qiskit.aqua.algorithms import VQE\n",
    "from qiskit.aqua.algorithms.adaptive.qaoa.var_form import QAOAVarForm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## High temperature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create an example system of two qubits that needs two extra qubits for purification. In this first example, we set $T=\\infty$, which corresponds to an inverse temperature $\\beta=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:11:13.932601Z",
     "start_time": "2018-11-19T20:11:13.928400Z"
    }
   },
   "outputs": [],
   "source": [
    "n_qubits = 2\n",
    "n_system = n_qubits * 2\n",
    "β = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define a simple Ising model with a weight matrix and set $p=1$ in QAOA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:11:13.948745Z",
     "start_time": "2018-11-19T20:11:13.937892Z"
    }
   },
   "outputs": [],
   "source": [
    "weights = np.array([[0,1],[0,0]])\n",
    "p = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on these parameters, we define the Ising Hamiltonian $H=\\sum W_{ij} \\sigma_i\\sigma_j$ (for the weight matrix defined above, $H=\\sigma_1 \\sigma_2$, whose minimum is reached when $\\sigma_1 \\neq \\sigma_2$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:11:13.984797Z",
     "start_time": "2018-11-19T20:11:13.950693Z"
    }
   },
   "outputs": [],
   "source": [
    "def pauli_z(qubit, coeff):\n",
    "    eye = np.eye((n_system))\n",
    "    return Operator([[coeff, Pauli(eye[qubit], np.zeros(n_system))]])\n",
    "\n",
    "def product_pauli_z(q1, q2, coeff):\n",
    "    eye = np.eye((n_system))\n",
    "    return Operator([[coeff, Pauli(eye[q1], np.zeros(n_system)) * Pauli(eye[q2], np.zeros(n_system))]])\n",
    "\n",
    "def ising_hamiltonian(weights):\n",
    "    H = reduce(lambda x,y:x+y,\n",
    "            [product_pauli_z(i,j, -weights[i,j])\n",
    "             for (i,j) in itertools.product(range(n_qubits), range(n_qubits))])\n",
    "    H.to_matrix()\n",
    "    return H\n",
    "\n",
    "Hc = ising_hamiltonian(weights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We prepare the intial state $|\\psi_0 \\rangle = \\sqrt{2 cosh(1/T)} \\sum_{z \\in {1, -1}} e^{- z/T} | z \\rangle_S \\otimes | z \\rangle_E$, with $E$ a temporary space used for purification purpose. It can be shown that tracing out this state over $E$ reproduces the state $\\rho \\propto e^{-H_m/T} $. We initialize the circuit first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:11:13.996055Z",
     "start_time": "2018-11-19T20:11:13.988499Z"
    }
   },
   "outputs": [],
   "source": [
    "qr = QuantumRegister(n_qubits * 2)\n",
    "cr = ClassicalRegister(n_qubits)\n",
    "backend = BasicAer.get_backend('qasm_simulator')\n",
    "circuit_init = QuantumCircuit(qr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we prepare the state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:11:14.017731Z",
     "start_time": "2018-11-19T20:11:14.010949Z"
    }
   },
   "outputs": [],
   "source": [
    "α = 2 * np.arctan(np.exp(- β/2))\n",
    "for i in range(n_qubits):\n",
    "    circuit_init.rx(α, qr[n_qubits+i])\n",
    "    circuit_init.cx(qr[n_qubits+i], qr[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The class `QAOAVarForm` takes an initial state as argument, which should be an object containing a method `construct_circuit` that builds and returns the initial state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:11:14.032266Z",
     "start_time": "2018-11-19T20:11:14.019967Z"
    }
   },
   "outputs": [],
   "source": [
    "class InitialState:\n",
    "    def __init__(self, β, n_qubits):\n",
    "        self.β = β\n",
    "        self.n_qubits = n_qubits\n",
    "    \n",
    "    def construct_circuit(self, mode, qr):\n",
    "        circuit_init = QuantumCircuit(qr)\n",
    "\n",
    "        alpha = 2 * np.arctan(np.exp(- self.β / 2))\n",
    "        for i in range(n_qubits):\n",
    "            circuit_init.rx(alpha, qr[n_qubits+i])\n",
    "            circuit_init.cx(qr[n_qubits+i], qr[i])\n",
    "            \n",
    "        return circuit_init"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to overload the original QAOA class to be able pass it an initial state as a parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:11:14.046628Z",
     "start_time": "2018-11-19T20:11:14.034301Z"
    }
   },
   "outputs": [],
   "source": [
    "class MyQAOA(VQE):\n",
    "    def __init__(self, operator, optimizer, init_state, p=1, operator_mode='matrix', initial_point=None,\n",
    "                 batch_mode=False, aux_operators=None):\n",
    "            self.validate(locals())\n",
    "            var_form = QAOAVarForm(operator, p, init_state)\n",
    "            super().__init__(operator, var_form, optimizer,\n",
    "                             operator_mode=operator_mode, initial_point=initial_point)\n",
    "        \n",
    "class MyQAOAVarForm(QAOAVarForm):\n",
    "    def construct_circuit(self, angles, quantum_register, classical_register):\n",
    "        if not len(angles) == self.num_parameters:\n",
    "            raise ValueError('Incorrect number of angles: expecting {}, but {} given.'.format(\n",
    "                self.num_parameters, len(angles)\n",
    "            ))\n",
    "        q = quantum_register\n",
    "        circuit = QuantumCircuit(quantum_register, classical_register)\n",
    "        if self._initial_state:\n",
    "            circuit += self._initial_state.construct_circuit('circuit', q)\n",
    "        else:\n",
    "            circuit.u2(0, np.pi, q)\n",
    "        for idx in range(self._p):\n",
    "            β, γ = angles[idx], angles[idx + self._p]\n",
    "            circuit += self._cost_operator.evolve(None, γ, 'circuit', 1, quantum_registers=q)\n",
    "            circuit += self._mixer_operator.evolve(None, β, 'circuit', 1, quantum_registers=q)\n",
    "        return circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We prepare the initial state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:11:14.069066Z",
     "start_time": "2018-11-19T20:11:14.049530Z"
    }
   },
   "outputs": [],
   "source": [
    "initial_state = InitialState(β, n_qubits)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We run the protocol to get the thermal state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:11:14.083442Z",
     "start_time": "2018-11-19T20:11:14.071127Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Begin QAOA...\n"
     ]
    },
    {
     "ename": "TypeError",
     "evalue": "construct_circuit() missing 1 required positional argument: 'qr'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-14-e34ea5f94119>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 13\u001b[0;31m \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_thermal_state\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mweights\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-14-e34ea5f94119>\u001b[0m in \u001b[0;36mget_thermal_state\u001b[0;34m(weights, p)\u001b[0m\n\u001b[1;32m      7\u001b[0m     \u001b[0mbackend\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mAer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_backend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'statevector_simulator'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m     \u001b[0mquantum_instance\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mQuantumInstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbackend\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mshots\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m     \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mqaoa\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mquantum_instance\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     10\u001b[0m     \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Results of QAOA\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/qiskit/aqua/qiskit/aqua/algorithms/quantum_algorithm.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, quantum_instance, **kwargs)\u001b[0m\n\u001b[1;32m     65\u001b[0m                 \u001b[0mquantum_instance\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_config\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     66\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_quantum_instance\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mquantum_instance\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 67\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_run\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     68\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     69\u001b[0m     \u001b[0;34m@\u001b[0m\u001b[0mabstractmethod\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/qiskit/aqua/qiskit/aqua/algorithms/adaptive/vqe/vqe.py\u001b[0m in \u001b[0;36m_run\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    283\u001b[0m                                       \u001b[0mvar_form\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvar_form\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    284\u001b[0m                                       \u001b[0mcost_fn\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_energy_evaluation\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 285\u001b[0;31m                                       optimizer=self.optimizer)\n\u001b[0m\u001b[1;32m    286\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    287\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_ret\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'num_optimizer_evals'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_eval_count\u001b[0m \u001b[0;34m>=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_ret\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'num_optimizer_evals'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/qiskit/aqua/qiskit/aqua/algorithms/adaptive/vq_algorithm.py\u001b[0m in \u001b[0;36mfind_minimum\u001b[0;34m(self, initial_point, var_form, cost_fn, optimizer, gradient_fn)\u001b[0m\n\u001b[1;32m    119\u001b[0m                                                                       \u001b[0mvariable_bounds\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mbounds\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    120\u001b[0m                                                                       \u001b[0minitial_point\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minitial_point\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 121\u001b[0;31m                                                                       gradient_function=gradient_fn)\n\u001b[0m\u001b[1;32m    122\u001b[0m         \u001b[0meval_time\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mstart\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    123\u001b[0m         \u001b[0mret\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/qiskit/aqua/qiskit/aqua/components/optimizers/cobyla.py\u001b[0m in \u001b[0;36moptimize\u001b[0;34m(self, num_vars, objective_function, gradient_function, variable_bounds, initial_point)\u001b[0m\n\u001b[1;32m     92\u001b[0m         \u001b[0msuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moptimize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnum_vars\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mobjective_function\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgradient_function\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvariable_bounds\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minitial_point\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     93\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 94\u001b[0;31m         \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mminimize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobjective_function\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minitial_point\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtol\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_tol\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"COBYLA\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moptions\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_options\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     95\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnfev\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.virtualenvs/aqua/lib/python3.7/site-packages/scipy/optimize/_minimize.py\u001b[0m in \u001b[0;36mminimize\u001b[0;34m(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)\u001b[0m\n\u001b[1;32m    604\u001b[0m                              **options)\n\u001b[1;32m    605\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'cobyla'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 606\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0m_minimize_cobyla\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconstraints\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    607\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'slsqp'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    608\u001b[0m         return _minimize_slsqp(fun, x0, args, jac, bounds,\n",
      "\u001b[0;32m~/.virtualenvs/aqua/lib/python3.7/site-packages/scipy/optimize/cobyla.py\u001b[0m in \u001b[0;36m_minimize_cobyla\u001b[0;34m(fun, x0, args, constraints, rhobeg, tol, maxiter, disp, catol, **unknown_options)\u001b[0m\n\u001b[1;32m    250\u001b[0m     xopt, info = _cobyla.minimize(calcfc, m=m, x=np.copy(x0), rhobeg=rhobeg,\n\u001b[1;32m    251\u001b[0m                                   \u001b[0mrhoend\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mrhoend\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0miprint\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0miprint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmaxfun\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmaxfun\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 252\u001b[0;31m                                   dinfo=info)\n\u001b[0m\u001b[1;32m    253\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    254\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0minfo\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0mcatol\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.virtualenvs/aqua/lib/python3.7/site-packages/scipy/optimize/cobyla.py\u001b[0m in \u001b[0;36mcalcfc\u001b[0;34m(x, con)\u001b[0m\n\u001b[1;32m    240\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    241\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mcalcfc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcon\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 242\u001b[0;31m         \u001b[0mf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    243\u001b[0m         \u001b[0mi\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    244\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0msize\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mizip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcons_lengths\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconstraints\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/qiskit/aqua/qiskit/aqua/algorithms/adaptive/vqe/vqe.py\u001b[0m in \u001b[0;36m_energy_evaluation\u001b[0;34m(self, parameters)\u001b[0m\n\u001b[1;32m    317\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0midx\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mparameter_sets\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    318\u001b[0m             \u001b[0mparameter\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mparameter_sets\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0midx\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 319\u001b[0;31m             \u001b[0mcircuit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconstruct_circuit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mparameter\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_quantum_instance\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbackend\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_use_simulator_operator_mode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    320\u001b[0m             \u001b[0mcircuits\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcircuit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    321\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/qiskit/aqua/qiskit/aqua/algorithms/adaptive/vqe/vqe.py\u001b[0m in \u001b[0;36mconstruct_circuit\u001b[0;34m(self, parameter, backend, use_simulator_operator_mode)\u001b[0m\n\u001b[1;32m    201\u001b[0m             \u001b[0;34m[\u001b[0m\u001b[0mQuantumCircuit\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mgenerated\u001b[0m \u001b[0mcircuits\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mHamiltonian\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    202\u001b[0m         \"\"\"\n\u001b[0;32m--> 203\u001b[0;31m         \u001b[0minput_circuit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_var_form\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconstruct_circuit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mparameter\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    204\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mbackend\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    205\u001b[0m             \u001b[0mwarning_msg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"Circuits used in VQE depends on the backend type, \"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/qiskit/aqua/qiskit/aqua/algorithms/adaptive/qaoa/var_form.py\u001b[0m in \u001b[0;36mconstruct_circuit\u001b[0;34m(self, angles)\u001b[0m\n\u001b[1;32m     58\u001b[0m         \u001b[0mcircuit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mQuantumCircuit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     59\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_initial_state\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 60\u001b[0;31m             \u001b[0mcircuit\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_initial_state\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconstruct_circuit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'circuit'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     61\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcircuit\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mqregs\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     62\u001b[0m             \u001b[0mq\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mQuantumRegister\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cost_operator\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnum_qubits\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'q'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: construct_circuit() missing 1 required positional argument: 'qr'"
     ]
    }
   ],
   "source": [
    "def get_thermal_state(weights, p):\n",
    "    Hc = ising_hamiltonian(weights)\n",
    "    print(\"Begin QAOA...\")\n",
    "    \n",
    "    optimizer = COBYLA()\n",
    "    qaoa = MyQAOA(Hc, optimizer, initial_state, p, operator_mode=\"matrix\")\n",
    "    backend = Aer.get_backend('statevector_simulator')\n",
    "    quantum_instance = QuantumInstance(backend, shots=100)\n",
    "    result = qaoa.run(quantum_instance)\n",
    "    print(\"Results of QAOA\", result)\n",
    "    \n",
    "    return result\n",
    "result = get_thermal_state(weights, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we reformat the final results, measure out the result, and plot the energy distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:11:44.601197Z",
     "start_time": "2018-11-19T20:11:14.085143Z"
    }
   },
   "outputs": [],
   "source": [
    "var_form = MyQAOAVarForm(Hc, p, initial_state)\n",
    "thermal_state = var_form.construct_circuit(result['opt_params'], qr, cr)\n",
    "for i in range(n_qubits):\n",
    "    thermal_state.measure(qr[i], cr[i])\n",
    "    job = execute(thermal_state, backend, shots=2000)\n",
    "results = job.result().get_counts(thermal_state)\n",
    "\n",
    "def get_energy(spin_configuration):\n",
    "    x = spin_configuration.reshape(-1, 1)\n",
    "    return np.sum([[-weights[i,j] * x[i] * x[j] for j in range(n_qubits)] for i in range(n_qubits)])\n",
    "\n",
    "list_spin_configs = np.array(np.concatenate([[list(spin_config)] * results[spin_config] for spin_config in results]), dtype=\"int\")\n",
    "list_spin_configs[list_spin_configs == 0] = -1\n",
    "list_energy = np.array([get_energy(spin_config) for spin_config in list_spin_configs])\n",
    "hist = plt.hist(list_energy, density=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two eigenvalues, i.e. possible energies, of our Hamiltonian $H=\\sigma_1 \\sigma_2$ are $E=-1$ and $E=1$. At infinite temperature ($\\beta=0$), they should be assigned an equal probability, which is the case in the histogram above. Let's repeat the experiment at a lower temperature."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Low temperature\n",
    "\n",
    "We set the inverse temperature to a high value. With this, we should get the lowest energy with a high probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "β = 5\n",
    "initial_state = InitialState(β, n_qubits)\n",
    "result = get_thermal_state(weights, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var_form = MyQAOAVarForm(Hc, p, initial_state)\n",
    "thermal_state = var_form.construct_circuit(result['opt_params'], qr, cr)\n",
    "for i in range(n_qubits):\n",
    "    thermal_state.measure(qr[i], cr[i])\n",
    "    job = execute(thermal_state, backend, shots=2000)\n",
    "results = job.result().get_counts(thermal_state)\n",
    "\n",
    "list_spin_configs = np.array(np.concatenate([[list(spin_config)] * results[spin_config] for spin_config in results]), dtype=\"int\")\n",
    "list_spin_configs[list_spin_configs == 0] = -1\n",
    "list_energy = np.array([get_energy(spin_config) for spin_config in list_spin_configs])\n",
    "plt.hist(list_energy, density=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The minimum energy eigenstate should now have a much higher probability. Try to repeat the experiment with different $\\beta$ to see the effect of the temperature on the prepared thermal state."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# References\n",
    "\n",
    "[1] Verdon, G., Broughton, M., Biamonte, J. (2017) [A quantum algorithm to train neural networks using low-depth circuits](https://arxiv.org/abs/1712.05304). *arXiv:1712.05304*. <a id='1'></a>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
